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Summary. The paper presents results for deriving closed-form analytic so- 
lutions of the non-relativistic linear perturbation equations, which govern the 
evolution of inhomogeneities in a homogeneous spatially flat multicomponent 
cosmological model. Mathematical methods to derive computable forms of the 
perturbations are outlined. 
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1. Introduction 

This paper deals with solutions of a system of differential equations describ- 
ing the evolution of density perturbations in the Universe. Before being more 
specific about the differential equations, we look at previous work setting the 
physical environment for deriving such equations. 

To formulate a quantitative picture of the evolution of perturbations outside 
the horizon it is necessary to solve the perturbed Einstein equations. The density 
contrast S(x) is not a gauge invariant quantity. For sub- horizon-sized pertur- 
bations, aX = Xphys <C if -1 , this fact is of little consequence, as a Newtonian 
analysis is sufficient; a denotes the scale factor and H the Hubble parameter, 
respectively. For super- horizon-sized perturbations, \ p h ys S> H -1 , the gauge 
transformation for 6(x) must be taken into account. In this case a full general 
relativistic treatment is required. So far two approaches are available for finally 
deriving the differential equation governing the growth or decay of gravitational 
instabilities in an expanding Universe. The metric perturbation approach was 
invented by Lifshitz (1946) with his derivation of the synchronous equations. 
This approach was later used by Bardeen (1980) to derive the comoving equa- 
tions. Hawking (1966) initiated an alternative approach to follow the evolution 
of the density perturbations in the expanding Universe by employing the gen- 
eral relativistic fluid flow equations. This method was used by Olson (1976) 
to derive equations which are equivalent to the synchronous equations, and the 
same method was further developed by Lyth and Mukherjee (1988) to provide 
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the comoving equations. 

In order to split up the energy density at a given point in spacetime into 
an average value plus a perturbation, one must specify a spacelike hypersurface 
on which the averaging is to be performed. The choice of coordinates in or- 
der to define the energy perturbation throughout spacetime is referred to as a 
choice of gauge (Bardeen 1980). Bardeen proposed to use the comoving gauge 
to deal with the general relativistic equations. Another choice is to use the 
synchronous gauge (Peebles 1980) which, however, is beset by the problem of 
introducing arbitrary "extra gauge mode" solutions. Actually there is an infinity 
of synchronous gauges. 

A choice of Lifshitz's or Hawking's approach and synchronous or comoving 
gauge having been made sets the framework for deriving differential equations 
describing the time evolution of the energy density perturbation. Here again 
one has the choice to subscribe to the perfect fluid approximation or the kinetic 
description (Peebles 1980). In the perfect fluid approximation, which is more 
physically transparent, the time evolution of the energy density perturbation 
is specified by a differential equation involving the pressure perturbation as an 
additional unknown. The equation is of second order in the comoving gauge, 
assuming that the pressure perturbation is known, and it has two solutions 
revealing a growing and decaying mode (Peebles 1980). Since the equation 
governing the amplitude 5(t) is of second order, there are two solutions. A given 
perturbation is expressed as a linear combination of 6+(t) and S-(t). At late 
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times, only the projection onto the growing mode may be important. Physically, 
the decaying mode corresponds to a perturbation with initial overdensity and 
velocity arranged so that the initial velocity perturbation eventually "undoes" 
the density perturbation. In synchronous gauge the equation is of third or higher 
order corresponding to the extra freedom of choosing the respective gauge. 

Many efforts have been made to derive approximate solutions of the equa- 
tions of rclativistic and non-relativistic perturbation theory, usually valid for 
early time large-scale and for late time small- scale inhomogeneities, respec- 
tively. The primary purpose of this paper is to derive closed-form analytic 
solutions of the non- rclativistic linear perturbation equations, which govern 
the evolution of inhomogeneities in homogeneous spatially flat multicomponent 
cosmological models. The mathematical method employed here can also be used 
to tackle the relativistic version of the equations involved, which will be done 
in a subsequent paper. These closed-form solutions are valid for irregularities 
on scales smaller than the horizon. They may be used for analytical interpo- 
lation between known expressions for short and long wavelength perturbations. 
In Section 2 we present the physical parameters relevant to the problem and 
the fundamental differential equations for arbitrary polytropic index 7, and a 
radiation or matter dominated Universe, respectively. Solutions of this funda- 
mental equation are provided in all subsequent Sections 3,4, and 5 for relevant 
parameter sets including the polytropic index ji and the expansion law index 77. 
To catalogue these solutions in closed-form, in all sections the theory of Meijer's 
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G- function will be employed. The following results are a continuation of the 
ones by Haubold, Mathai, and Mucckct (1991), Mathai (1989), and Nurgaliev 
(1986). 

2. Some Parameters of Physical Significance 

The growth of density inhomogeneities can begin as soon as the Universe is 
matter-dominated. This would be also the case in a Universe dominated by non- 
interacting relic Wimps. Baryonic inhomogeneities cannot begin to grow until 
after decoupling because until then, baryons are tightly coupled to the pho- 
tons. After decoupling, when the Universe is matter-dominated and baryons 
are free of the pressure support provided by photons, density inhomogeneities 
in the baryons and any other matter components can grow. Actually the time of 
matter-radiation equality is the initial epoch for structure formation. In order to 
fill in the details of structure formation one needs " initial data" for that epoch. 
The initial data required include, among others, (1) the total amount of non- 
relativistic matter in the Universe, quantified by Oo, and (2) the composition of 
the Universe, as quantified by the fraction of critical density, 0, = pi/p c , con- 
tributed by various components of primordial density perturbations (i= baryons, 
Wimps, relativistic particles, etc.). Here the critical density p c is the total mat- 
ter density of the Einstein-de Sitter Universe. Speculations about the earliest 
history of the Universe have provided hints as to the appropriate initial data: 
flo = 1 from inflation; 0.015 < fis < 0.15 and SIwimp ~ 0.9 from inflation, 
primordial nucleosynthesis, and dynamical arguments. In the following we as- 
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sume that ^2^i — 1,0, = const.. The cosmological medium is considered 
to be unbounded and that there may be a uniform background of relativistic 
matter. It is further assumed that the matter is only slightly perturbed from 
the background cosmological model. This assumption may give a good descrip- 
tion of the behavior of matter on large scales even when there may be strongly 
nonlinear clustering on small scales (primordial objects). Also assumed is that 
matter can be approximated as an ideal fluid with pressure a function of density 
alone. It consists of i components having the densities pi and velocities of sound 
01 = dPi/dpi oc p 7i_1 , when an equation of state 

P l oc p^ (1) 

has been taken into account. The i components of the medium are interrelated 
through Newton's field equation A</> = 4ttG Pi, containing the combined den- 
sity Pi °f au components. Superimposed upon an initially homogeneous and 
stationary mass distribution shall be a small perturbation, represented by a sum 
of plane waves 

5 = 5 -£i = 5i{t)e ikx , (2) 
P 

where k = 2ira/\ defines the wave number k and A is the proper wavelength. In 
the linear approach the system of second order differential equations describing 
the evolution of the perturbation in the non-relativistic component i is 

j=i 
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where for all perturbations the same wave number k is used. This equation is 
valid for all sub-horizon-sized perturbations in any non-relativistic species, so 
long as the usual Friedmann equation for the expansion rate is used. Although 
the following investigation of closed-form solutions of the fundamental differen- 
tial equation governing the evolution of inhomogcneities in a multicomponent 
cosmological medium is quite general, it will be presented within the context of 
the inflationary scenario. Therefore in the following we consider an Einstein-de 
Sitter Universe with zero cosmological constant. After inflation the Friedmann 
equation for the cosmological evolution reduces to 



The continuity equation gives p oc a~ 3 and from cq. (4) one has a oc i 2 / 3 
and pi = Qi/(6nGt 2 ). During the expansion the Hubble parameter changes as 
H = rjt^ 1 , where the expansion law index is rj = ^ in the radiation-dominated 
epoch and 77 = | in the matter-dominated epoch respectively. The wave number 
is proportional to a -1 so that k 2 (3f = kft 21 ^ 1 ' 71 " 1 ^ in cq. (3), where now the 
constants fcj come from both the wave vector and the velocity of sound. If k{ = 
the sound velocity fa equals zero, in which case the adiabatic index 7$ loses 
its sense. Defining the parameter on = 2(2 — 77 — 7,) that absorbs the adiabatic 
index as well as the type of expansion law we can write for eq. (3) using eq. 



(4): 
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and dots denote derivatives with respect to time. We introduce a time oper- 
ator A = td/dt and change the dependent variable, the density perturbation 
Si(t), dealing with the function <j>i instead of Si by setting <5,(t) = t a <&i(t). The 
equation (5) for the function $j is then given by 

2 m 
3=1 



where 



= r a S u h = k 2 t a > - a 2 , on = 2(2 - n - 7i), ol = - (^j^) ■ (6) 



Observing that — 1 and operating on both sides of the differential equation 
for the $i by A 2 we obtain the fundamental equation 



A 4 $ l + A 2 (b^)-l(A 2 ^ + b l ^) = -^jrb j n 3 <i> J1 i = l,...m, (7) 

(Mathai 1989). As indicated above for r] the values | and \ are significant and 
in what follows it will be assumed that 2 > 7$ > |. We have chosen the range of 
values of 7, for physical as well as mathematical reasons as will be evident later 
in the analysis of eq. (7). For some of these parameter values we will consider 
a multicomponent medium. Consider the special case 

&1 = ^-fcD! 



, =6 = fcV _(^-P 2 
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of the fundamental equation (7). This is the case considered in Haubold, Mathai, 
and Muecket (1991) in detail. In the present discussion we will use the same 



notations as in Haubold, Mathai, and Mueckct (1991). Let k = in (8). This is 
case 4.1 of Haubold, Mathai, and Mueckct (1991). For this case the parameters 
are the following: 



K,b* = ± 



f (2*7-1 
\ 4a? 



,2 1 2 



at 



(2r?-l) 2 2 
4 3 



a 1; a 2 



-1± 



(2*7 -I) 2 , 2 2 
^^ + 3~3^ 



Table 2.1 gives these parameters fort)= | ,| and 7, = §, 1, |, 1, 2. 



(9) 
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Table 2.1 



V, 7* 




a 


h* h* 
°1> °2 


h* h* 
°3> °4 


a*,a2 /or fii = ± 


2 4 

3 3 





l 

6 








- 1 

3 L 


2 
3 


1 
6 


±3 


±1 


_1±^3 


2 2 

3 3 


4 
3 


1 
6 




±1 


/TT7 

1 ^ 8 


2 5 

3 3 


2 
3 


1 
6 


±4 


±1 




- 2 

3 Z 


4 
3 


1 
6 


±1 


±^ 

8 


/To 

-1± vM 

8 


1 4 

2 3 


1 
3 








±V6 


-l±v/3 


|1 


1 








V 3 


-^73 


1 2 

2 3 


5 
3 








_|_ 

5 


-l±f 


1 5 

2 3 


1 

3 








±V6 


-1± V3 


± 2 

2 Z 


-1 








V 3 





When evaluating the solutions explicitly one needs the parameter differences. 
These are given in Table 2.2. 
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Tabic 2.2 



V 


7i 


°1 


b 2 


b 3 


o 4 


b l _ b 2 


°1 _ °3 


°1 _ °4 


°2 _ °3 


°2 — °4 


3 - 4 


2 
3 


1 


1 

4 


1 

4 


5 

4 


5 
4 


1 

2 


— 1 


3 
2 


3 
2 


1 


5 
2 


2 
3 


2 
3 


1 

8 


1 

8 


5 

8 


5 
8 


1 

4 


1 

2 


3 
4 


3 
4 


1 
2 


5 
4 


2 
3 


5 
3 


4 


I 

~4 


5 
4 


5 

— 4 


I 

2 


-1 


3 
2 


3 

— 2 


1 


5 
2 


2 
3 





1 

8 


1 

8 


5 

8 


5 

8 


1 

4 


1 

2 


3 
4 


3 
4 


1 
2 


5 
4 


1 
2 


4 
3 








V6 


-V6 





-V6 


V6 




V6 


2V6 


1 
2 


1 








V 3 


V 3 





[2 

V 3 


[2 

V 3 


[2 

V 3 


[2 

V 3 


Z V 3 


1 
2 


2 
3 








5 


5 





5 


5 


5 


Ve 

5 


9 76 
z 5 


1 
2 


5 
3 










-VS 





-V6 








2V6 


1 
2 


2 








V 3 


~V 3 





-V 3 


V 3 


-V 3 


V 3 


2 V 3 
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3. Solution when oi\ = 

Note that when r\ = | and 7i = § one can have on = in equation (6). The 
exact solution given in equation (5.3) of Haubold, Mathai, and Mucckct (1991) 
is for the situation on ^ 0. Hence this case needs separate discussion. When 
ai = the fundamental equation (7) reduces to the form 

'(2r? - l) 4 2 (2?7 - l) 2 



A 4 $x 



(2TJ-11 2 f 
2 3 1 



2„ (2r]-l) 2 2 
l : — 
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$i = 0. 



(10) 



In this case the general solution is of the form 



$i = cit dl + c 2 t a2 + c 3 i 03 + 



<^3 I 



(11) 



where c\, C2, C3, C4 are arbitrary constants and di, cfo, cfa, c?4 are the roots of the 
equation 

'(2r,-l) 2 , 2 



+ 



(2r?-l) 4 2 (2r?-l) 2 [2 (2^-l) 2 2. 
T^ + 3 4 + l3 Sil 4 3^ 



0.(12) 



They can be seen to be the following: 

d\ , d,2 = ± 



(*L^ + I_M + J ( i + *?)_H ni 

4 3 2 1 v 3 2 7 3 



d 3 ,d 4 = ± 



(2V-11 i__kj _r 1 ^ 2 

4 3 2 P3 2 ; 3 



(13) 



(14) 
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4. Solution when the parameters differ by 
integers 

The general solution given in Haubold, Mathai, and Muecket (1991), equa- 
tions (5.2) and (5.3), are for finite values of t and for the cases that the 6*'s 
do not differ by integers. From Table 2.2 note that at two points, namely 
(r) = §,7, = §) and (77 = §,7j = 1) one has b\ - = -1 and b\ - b 2 = -1. 
This means that G\ and G4 of (5.3) of Haubold, Mathai, and Muecket (1991) 
need modifications. At all other points the G/s remain the same. 

4.1. Modifications of the G'jS for the cases (rj = |,7j = 1) 

and (77 = |,7i = §) 

Consider, for x = k\t ai /a\, ot\ 7^ 0, 

1 r r(l + S )T(-al - s)T(-a* 2 - S ) X -°d S 

z™Jl r(|- s )r(-i- s )r(|- s ) ' ( > 

Note that a zero coming from r(— \ — s) at s = —j coincides with a pole coming 
from T(i + s) at s — —\. This can be removed by rewriting as follows: 

G i = 7T— / ^75 777^ 777779 ^ ■ ( 16 ) 



2« J l r(| - s )T(i - S )r(| - a) 

Evaluating (16) as the sum of the residues at the poles of T(| + s) one has 



Gi = -x 



r(- a * + |)r(-^ + |) 
r(f)r(2)r(M } 



< * 5 * 5 10 14 

x 2 F 3 (-a* 1 + -,-a* 2 + -;—,2,—;-x). (17) 
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Evaluating G4 also the same way one has the following: 



r _ r l,2 ( ,a*+l,a*+l\ 

J_ /■ r(-f + s )r(- a |- s )r(-^- a ) _ 
2«A r(f- s )r(f- s )r(-i- s ) x rfs 

- 2™ A r(l- s )r(f- s )r(-i- s ) x ds 
_ 4 r(-i-a;)r(4-a^) 
r(|)r(2)r(-i) 

x iFz{-\-a\,-\-al-, l -,2- l ---x). (18) 

Thus the complete solution in this case is of the form 

$1 = c 1 G 1 +c 2 G 2 + c 3 G 3 + c i G i , (19) 

where c\, c 2l c 3l C4 are arbitrary constants, G\ and G4 are given in (17) and (18) 
respectively and G 2 and G 3 are given in equation (5.3) of Haubold, Mathai, and 
Muecket (1991). 

5. Solution near oo for the parameter values of 
table 2.1 

Except for the points (77, 7i) = (|, |), (|, 2) in all other cases some of the 
parameters differ by integers. In five cases of Table 2.2 the poles in the integrand 
can be up to order 2 and in two cases the poles can be of order up to 3. The 
general solution near 00 is of the form 

$1 - fiF! + f 2 F 2 + f 3 F 3 + f 4 F 4 , (20) 
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where /i, f 2 , /3, /4 are arbitrary constants and Fj are given in the equation (5.6) 
of Haubold, Mathai, and Muecket (1991). For example, 

Fi = G 4 2 i (x ) ,x = fc?t°Va?,ai + 0. (21) 

F2 is available from F\ by interchanging a\ and a^. -F3 and F4 are available 
from Fx by replacing x by xe l7r and se" 1 " respectively. We will evaluate F\ for 
the parameter values which differ by integers. 

5.1. F 1 for (r/=§,7; = §) 
In this case one has 

j_ /■ rwrwye + ^-Ve + ^-ai-,) 

2tt*./ l r(l + a* + s) v ; 

Note that at s = 0, —1, —2, ... the integrand has poles of order 2 each. All other 
poles arc of order 1 each. Thus 

F 1 =H 1 +H 2 + H 3 , (23) 

where H 2 and H 3 are the sums of the residues at the poles of r(-\/6 + s) and 
r(— \/6 + s) respectively and Hi is the sum of the residues at the poles of T 2 (s). 
For all the cases of simple poles including H 2 and H3 one has the result as 
follows: 

F(-at + b*) 



Rj = X h 'i 



KUm - b*) 



T(l + a*-b*) 

x 2 F s {-a\ + b*,-a* 2 + b*; 1 - b\ + b*, 1 - b* 4 + b*; -x)(24) 

where II' denotes the absence of the gamma F(&£ — b* k ) and # denotes the 
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absence of the element 1 — b* — b*. Thus one has, 

x 2^3 (~at + V6, -a* 2 + V6; 1 + VE, 1 + V6, 1 + 2^6; -a:) (25) 



and 



H 3 = X -^T\V6)T(2V6) ^-f:f^ 

T(l + a* 2 + V6) 

x 2-F3(-ai - V6, -al - 1 - V&, 1 - 1 - 2\/6; -a;). (26) 



Now consider the evaluation of Note that for all parameter combinations 

{V,li) = (|>|)' (I' 1 )' (I'D' (I'D' (I' 2 ) we have & i = = h *2- Hcnce we can 
write H\ in general terms as follows, 

Hi = j_r THs)m+sn b i + sn-at-s) 

2m J L T(l+a* 2 + s) V 1 

The general technique of evaluating such integrals is available in Mathai and 
Saxena (1973). The solution is the following 
1 



* = £ 



- (lux) {2i/>(v + 1) + V(&3 - f) + ^(62 - v) 



v=0 

+i>{-al + v)-i,{l + a* 2 -v)}] 

where the 63 and b\ are given in Table 2.1 for the specific parameter combina- 
tions and ip() is a psi function. Note that the first part in (27) can be written 
as a 2^3(-Oi) -o-2'i 1) 1 ~ b h 1 - &2; 
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5.2. Fx for ( r/ = |, 7i = l) and (?? = |,7^ = |) 



In these two cases Fi is of the following form: 

1 f T(\ + s)T(-\ + s)r(| + s)r(-f + a )r(-ai - S K s ds 



2niJ L {s + l) { s-l)T(l + a* 2 + s) 



2iri J L T(l + a* 2 + s) 

i r r 2 (| + s )r 2 (-i + s )r(- a t- s ) 



There are poles of order one each at s = — \ and s = | respectively and poles 
of order 2 each at s = — | — v,s = j — A, z/ = 0, 1, . . . , A = 0, 1, . . . respectively. 
Thus Fx is evaluated as the sum of the residues at all these poles. Let us denote 
these residues by Ri, i?2, R3 and R4 respectively. Then 

F 1 =R 1 +R 2 +R z + R Al (30) 

where 



i?i = residue at s = — — 

4 

r 2 (-i)r(-q; + i) , 2 r 2 H)r(- a * + 1) , 

(i)r(| + a5) 3 r(| + a5) x ' 

5 

i? 2 = residue at s = — 

2r 2 (f)r(-at + |) _ s 

x 4 ; 



InT -I- J.ihdj — I — 1 — I — y.ihi — 

2 ' rv 1 4 



5 f 3 
i? 3 = \ -Ikx + 2tJj(iv + 1) + 2^(-- - v) +tp(-al + - A +v) 



5,1 1 



^ (1 + a ^I^ ) + ^TT + (|T^) 
r 2 (-|-^)r(-at + f + v) x" 

X T{l + al-\- V )(l + V )(\ + v) H) 2 ' 

1 00 r 3 

i? 4 = ^ J _ln. T + 2ip(v+ 1) + 2i/>(- - v) 



18 



Y^l-u)Y(-al-\ + V ) x » 
(\-v)(-l-v)T{\ + a* 2 -v)W K ' 

As before F 2 is available from the F\ of (30) by interchanging a\ and a 2 ■ F 3 and 
F 4 are available from F\ by replacing x by e t7r x and e _ "x respectively. This 
completes the evaluation of $1 for all finite values as well as for values near oo 
for all the parameter combinations given in Tables 2.1 and 2.2. 

6. Conclusion 

We have presented closed-form solutions of the non-relativistic linear perturba- 
tion equations which govern the evolution of inhomogeneities in a spatially flat 
multicomponent cosmological medium. The general solutions are catalogued 
according to the polytropic index ji and the expansion law index rj of the mul- 
ticomponent medium . All general solutions are expressed in terms of Meijer's 
G-function and their Mellin-Barnes integral representation. The proper use of 
this function simplifies the derivation of solutions of the linear perturbation 
equations and opens ways for its numerical computation. In this regard the 
paper leaves room for further work. 
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